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ABSTRACT 


The  differentiation  matrix  for  a  Daubechies-based  wavelet  basis  defined  on  an  in¬ 
terval  will  be  constructed.  U  will  be  shown  that  the  differentiation  matrix  based  on 
the  currently  available  boundary  constructions  does  not  maintain  the  superconver¬ 
gence  encountered  under  periodic  boundary  conditions. 
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□  □ 


1  Introduction 


Super  convergence  is  a  property  of  wavelet  methods,  Daubechies,  see  [4],  and  spline, 
see  [10],  based,  when  periodic  boundary  conditions  are  assumed.  That  is,  for  Daubechies- 
based  methods  under  the  assumption  of  periodicity  it  was  proven  in  [6]  that  the 
differentiation  matrix  is  accurate  of  order  2M,  even  though  the  highest  degree  poly¬ 
nomial  which  can  be  reconstructed  exactly  in  the  wavelet  subspace  is  of  degree  M  —  1, 
where  M  is  the  number  of  vanishing  moments  of  the  wavelet.  For  spline-based  meth¬ 
ods,  again  assuming  periodicity,  it  was  proven  in  [7]  that  the  differentiation  matrix 
is  accurate  of  order  2n  -|-  2,  even  though  one  can  only  construct  splines  of  order  n 
exactly  in  the  underlying  subspace.  For  Daubechies  and  spline-based  wavelet  systems 
this,  roughly,  doubling  of  the  differentiation  accuracy  is  known  as  superconvergence, 
a  name  borrowed  from  the  finite  element  literature. 

When  periodicity  is  no  longer  assumed  and  one  defines  wavelets  on  an  interval 
one  of  the  goals  when  building  boundary  functions  should  be  to  maintain  the  su¬ 
perconvergence  of  the  differentiation  matrix  across  the  entire  interval.  It  has  been 
proven  by  Gottlieb,  et.  al.  [5],  that  the  superconvergence  encountered  in  finite  el¬ 
ement  methods  under  the  assumption  of  periodicity  can  not  be  maintained  on  the 
interval  if  there  are  characteristics  leaving  the  domain.  Furthermore,  it  was  shown  in 
[7]  that  for  spline-based  wavelet  systems  that  superconvergence  is  lost  at  the  bound¬ 
aries  when  the  boundary  functions  are  truncated  B-splines.  The  goal  of  this  paper 
is  to  explore  the  accuracy  of  the  differentiation  matrix  for  Daubechies-based  wavelet 
systems  defined  on  an  interval.  In  particular,  the  differentiation  accuracy  of  the  dif¬ 
ferentiation  matrices  constructed  from  the  currently  available  boundary  constructions 
for  Daubechies  wavelets  defined  on  an  interval  will  be  explored.  The  differentiation 
matrix  for  the  Daubechies  4  wavelet  basis  using  the  boundary  functions  defined  in 
[3]  will  be  calculated  explicitly  and  the  accuracy  will  be  found.  From  this  explicit 
example  one  can  see  a  necessary  condition  for  superconvergence  to  be  maintained 
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with  other  boundary  constructions. 


We  begin  by  defining  the  differentiation  matrix,  P,  as  the  product  of  three  ma¬ 
trices:  Tf  =  CDC.  The  three  matrices  are  defined  as  follows: 

•  The  first  matrix  C  is  a  quadrature  matrix  mapping  from  samples  in  the  physical 
space  to  approximate  scaling  function  coefficients  at  the  finest  scale:  s  =  Cf. 
C  can  be  derived  from  the  moments  of  the  scaling  function.  (7,  also,  can  simply 
be  the  inverse  of  the  third  matrix:  C  =  (<7)“’. 

•  The  second  matrix  D  maps  from  the  scaling  function  coefficients  of  a  function 
to  the  scaling  function  coefficients  of  the  derivative  of  the  same  function.  For 
convenience  I  will,  henceforth,  refer  to  this  matrix  as  the  ‘derivative  projection’ 
matrix. 

•  The  third  matrix  C  is  a  quadrature  matrix  mapping  from  approximate  scaling 
function  coefficients  to  approximate  point  values  in  the  physical  space:  /  =  Cs. 
C  can  be  derived  from  samples  of  the  scaling  function  or  can  simply  be  the 
inverse  of  the  first  matrix:  C  =  C~^. 

This  paper  contains  the  following  sections. 

1.  Introduction 

2.  Boundary  Functions: 

An  outline  of  the  boundary  functions  constructed  in  [3]  is  given. 

3.  Derivative  Projection  Matrix  D: 

The  entries  of  D  are  defined  and  calculated.  D  has  a  banded  structure  due 
to  the  local  support  of  the  basis  functions,  but  D  is  not  antisymmetric.  The 
number  of  nonzero  entries  in  each  row  of  D  bounds  the  maximum  differentiation 
accuracy  for  the  corresponding  row  in  the  differentiation  matrix.  Furthermore, 
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the  number  of  nonzeros  entries  in  the  first  few  and  last  few  rows  of  D  are  related 
to  the  overlap  of  the  boundary  functions. 

4.  Approximating  Scaling  Function  Coefficients: 

An  outline  of  scaling  function  coefficient  approximation  is  given. 

5.  Quadrature  Matrix  C  Using  Moments: 

C  performs  the  mapping,  C  :  f  s.  The  D4  wavelet  has  two  vanishing 
moments  which  yields  well-conditioned  2-point  quadrature  matrices  and  ill- 
conditioned  3-point  quadrature  matrices. 

6.  Quadrature  Matrix  C  Using  Samples: 

A  A  ^  A 

C  performs  the  mapping,  C  :  s  /.  For  the  D4  wavelet  the  matrix  C  has 
three  entries  in  each  row  and  is  banded.  C  is,  however,  ill-conditioned. 

7.  Examples  and  Order  of  Accuracy; 

The  matrix  D  is  fixed.  There  are,  however,  many  choices  for  the  quadrature 
matrix  producing  many  versions  of  2>.  In  all  cases  superconvergence  is  lost  at 
the  boundaries.  Furthermore,  choosing  among  the  quadrature  matrices  which 
not  ill-conditioned  one  gets  a  differentiation  matrix  which  is  essentially  full.  A 
banded  V  can  be  constructed  from  a  combination  of  the  two  types  of  quadrature 
matrices  C  and  C.  The  bandwidth  from  such  a  construction  is  large,  9,  for 
a  D4  wavelet  and,  again,  the  superconvergence  is  lost  at  the  boundaries.  In 
addition,  all  of  these  matrices  display  directional  dependency.  This  directional 
dependency  can  be  ‘fixed’  by  reversing  the  coordinate  system  and  averaging  the 
appropriate  differentiation  matrices. 

8.  Conclusion: 

The  differentiation  matrices  for  wavelets  on  an  interval  using  the  boundary  func¬ 
tions  defined  in  [3]  do  not  maintain  the  superconvergence  encountered  when 
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periodic  boundary  conditions  are  assumed.  By  simply  counting  the  overlap  of 
boundary  functions  in  other  boundary  constructions  it  can  be  seen  that  all  cur¬ 
rently  available  boundary  constructions  fail  to  maintain  the  superconvergence 
encountered  under  periodic  boundary  conditions. 
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2  Construction  of  Boundary  Functions 


This  section  is  a  brief  outline  of  the  boundary  function  construction  proposed  by 
Cohen,  Daubechies,  and  Vial,  see  [4]. 

The  goal  is  to  build  a  wavelet  basis  on  an  interval  where  the  scaling  functions  and 
wavelets  away  from  the  boundaries  are  the  usual  Daubechies  scaling  functions  and 
wavelets.  At  the  boundaries,  boundary  scaling  functions  are  constructed  such  that 
polynomials  of  order  up  to  the  number  of  vanishing  moments  of  the  wavelet  can  be 
reproduced  exactly  across  the  entire  interval. 

The  boundary  function  construction  begins  by  building  independent,  but  not  or¬ 
thogonal,  functions  as  follows: 

2N-2  /  \ 

^  j  <f>ix  +  n-  N  +  1),  (1) 

where  4>{x)  is  the  usual  Daubechies  scaling  function  and  N  is  the  number  of  van¬ 
ishing  moments  of  the  associated  wavelet.  Note  that  these  functions  are  compactly 
supported,  and  their  supports  are  staggered,  i.e., 

supp(<^*')  =  [0,2iV  —  1  —  fc].  (2) 

The  staggered  support  yields  independence,  and  the  boundary  functions  are  then 
defined  by  simply  orthonormalizing  these  functions  by  a  Gram-Schmidt  method. 

Coefficients  and  can  be  found,  see  [4],  such  that  the  boundary  functions  are 
defined  recursively  as  follows  beginning  with  the  left-hand  side  boundary  functions: 

N-\  N+2k 

XI  X  -  m),  (3) 

/=0  m=/V 

and  the  equation  defining  the  right-hand  side  boundary  functions  is, 

N-l  N+2k 

</>k(^)  =  V2  X  ^  X  <m<^(2x  +  m+l).  (4) 

1=0  m=N 

The  following  figure  is  a  plot  of  the  four  boundary  scaling  functions  for  the 


Daubechies  4  wavelet. 


LHS  Boundary  Function  0 


LHS  Boundary  Function  1 


Figure  1:  The  boundary  functions  for  the  Daubechies  4  scaling  function. 


3  The  Derivative  Projection  Matrix  D 


Let  us  recall  the  origin  of  the  matrix  D.  First,  we  begin  with  a  function  f{x)  € 
L^{R).  Next,  we  want  to  approximate,  or  discretize,  this  function  in  a  scaling  function 
subspace  Vq  C 

PvJ{x)  =  "'£sM^).  (5) 

k=0 

where  {6/t(x)}  denotes  the  bases  functions  of  this  finest  scale  discretization  and  d 
is  the  dimension  or  number  of  degrees-of-freedom  of  the  subspace.  In  the  case  of 
wavelets  with  periodic,  boundary  conditions  =  <^*(x),  the  usual  Daubechies 

scaling  functions,  for  all  k  £  {0,...,d—  1}.  In  the  current  scenario,  however,  of 
wavelets  on  an  interval  bk{x)  is  different  at  the  boundaries; 

bk{x)  =  <f>^{x) 


for  k  =  0, N  —  \, 

for  A'^  <  k  <  d  —  2N  +  1  and. 


bk{x)  =  (pk{x) 

bkix)  =  (f>^ix) 


for  k  =  d  —  2N  +  2,...d  —  1,  where  N  is  the  number  of  vanishing  moments  of  the 
wavelet.  Under  the  orthonormality  of  the  current  basis  set,  <  6^,  6/  >=  Sk,h  where. 


<g,h  >=  j^g{x)h{x)dx, 

and  fj  denotes  integration  over  the  interval,  we  get. 


(6) 


Sk  =<  bk,f  >  . 


(7) 


We  now  differentiate  the  approximation  Pv„f{x)  to  get, 


where  b{x)  =  ^6(a:).  The  derivative  takes  Pvof{x)  out  of  Vo-  Projecting  back  into  Vo 
we  get, 


Pvo^PvJix)  =  f2<  >  b,(x), 


1=0 


dx 


or, 


d-i  d-i 


Pvo^PvoKx)  = 


(9) 


(10) 


1=0  k=0 

The  elements  <  6*,  6/  >  comprise  the  derivative  projection  matrix  D.  If  we  let  .7 
denote  the  vector  of  the  scaling  function  coefficients  for  k  =  0,...,d  —  1,  then  1) 
maps  from  the  scaling  function  coefficients  of  a  futiction  at  the  finest  scale  to  (he 
scaling  function  coefficients  of  the  derivative  of  the  same  function: 


D  : 


:  s  —4  5. 


D,  in  effect,  contains  all  the  derivative  information  in  a  wavelet  expansion.  The 
elements  of  D  contain  the  numerical  values  of  the  derivative  of  each  scaling  function 
projected  onto  all  scaling  functions  in  the  basis.  Due  to  the  compact  support  of 
the  bases  functions,  the  matrix  D  has  a  band  diagonal  form.  The  following  matrix 
illustrates  the  form  of  D  for  the  wavelet  basis  with  boundary  functions  defined  as 
in  the  previous  section: 


Po,o 

P\,o 

®2,0 

0 

0 

0 

0 

0 

0 

0  ■ 

Po,\ 

®2,1 

®3.1 

0 

0 

0 

0 

0 

0 

^o;i 

^,".2 

ro 

r2 

0 

0 

0 

0 

0 

0 

^1(3 

r-\ 

^’o 

ri 

r2 

0 

0 

0 

0 

0 

0 

r-2 

r-i 

ro 

n 

r2 

0 

0 

0 

0 

0 

0 

r-2 

r-i 

ro 

rj 

r2 

0 

0 

0 

0 

0 

0 

r-2 

r-i 

ro 

ri 

0 

0 

0 

O' 

0 

0 

r-2 

r-i 

ro 

0 

0 

0 

0 

0 

0 

“4,1 

Pu 

PoA 

0 

0 

0 

0 

0 

0 

0 

“5,0 

P\,0 

Po,o  . 

where  the  coefficients  r  are  the  derivative  projections  due  to  scaling  function  inter 
action,  the  coefficients  p  are  the  projections  of  the  derivatives  of  boundary  functions 
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projected  onto  boundary  functions,  and  the  coefficients  a  and  (3  contain  the  interac¬ 
tions  between  derivatives  of  boundary  functions  and  scaling  functions. 

In  the  following  subsections  the  entries  of  D  will  be  found. 

3.1  The  Equations  for  the  Derivative  Projections 

To  find  the  elements  of  D  one  nee<ls  to  find  the  interaction  of  the  derivative  of  each 
basis  function  with  every  otluT  basis  function.  Due  to  the  local  support  of  the  ba.ses 
functions,  the  left-hand  side  boimdary  functions  and  the  right-hand  boundary  func¬ 
tions  due  not  interact.  (Assuming,  of  course,  that  one  does  not  choose  an  unreason¬ 
ably  small  domain  relative  to  the  support  of  the  scaling  functions.) 

Each  entry  in  the  matrix  D  must  be  found  separately,  but  the  calculations  are 
straightforward.  A  sampling  of  the  relevant  derivations  will  be  given  below. 


3.1.1  Derivation  of  the  LHS  Equations 

Recall  that  the  LHS  boundary  functions  4>k{x)  are  generated  by, 

N-\  V+2it 

<f>i{x)  =  V2^hl4,f{2x)  +  ^/2  ^  hj;^^^{2x-m).  (11) 

/=o  m=yv 

Differentiate  the  above  equation  to  get, 

j  N+2k 

— =  2V2  +  2v^  y;  kt,J(2x  -  m).  (12) 

1=0  m=N 

Multiply  <^jt(x)  by  (f>p{x)  and  integrate  to  get, 

N-\N-\  . 

/  =  4  J:  5:  ^l'{2x)<l>^{2x)dx  (13) 

t=0  t=0 

N+2kN-\  . 

+  4  ^  I  (2x)(^(2i  -  m)dx 

m=N  i=0 
N-^  N+2p 

^  E  /  ^j'{2x)d>{2x  -  q)dx 

1=0  q=N  ■'' 

N-i-2k  N+2p 

+  4  ^  f^k.m^p.q  /  <^(2x  -  m)4>i2x  -  q)dx. 

m=N  q=N  ^ 
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Rename  the  integrals  as, 


Pk,p  = 

J^^ki^)4>p{x)dT, 

(14) 

= 

j  ij>{x  —  m)<j)^{x)dx, 

(15) 

4l'{x)<f>{x  -  q)dx, 

(16) 

^m,q  ~ 

j  (j>{x  —  m)(f>{x^  —  q)dx, 

(17) 

to  get  the  following  system  of  equations: 

N-\N~l  N-lN+2k 

=  E  E  ^ +  E  E + 

/  t  t=0  tn=N 

N-iN+2p  N+2kN+2p 

E  E  WJt,+  E  E  '■L'-L’-.m-  (18) 

1=0  q=S  m  =  N  q=N 

The  numerical  values  of  the  coefficients  {r}  were  found  in  [2],  and  the  numerical 
values  of  the  coefficients  o^,  and  will  be  found  in  the  following  subsections. 

3.2  The  Coefficients  pf;,,  for  i  =  0, 1 

In  this  subsection  we  will  find  the  corner  diagonal  elements  of  D.  These  elements 
represent  the  projection  of  the  derivative  of  each  boundary  function  onto  itself. 

3,2.1  The  Coefficients  p^g  and  pf  , 

For  simplicity,  allow  the  lower  limit  to  be  0  and  the  upper  limit  be  oo  in  the  definition 

of  P^k- 

Pk,k  =  /  <f>k{x)(f>i{x)dx. 

Jo 

Apply  integration  by  parts  to  this  equation  to  get, 


Pi.1.  =  ((‘f(^))'IS°  -  / 
^0 


^Pk,k  — 
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(21) 


<j>k{x)  is  0  at  X  =  oo  leaving  only  the  value  at  the  lower  limit, 

Pk,k  -  2  • 

That  is,  pj^i^  is  determined  by  the  value  of  the  boundary  function  at  the  boundary. 
Before  evaluating  the  boundary  functions  at  the  boundary  it  will  be  shown  that  Pq  q 
and  p^i  are  related  by  a  constant  multiple. 

3.2.2  The  Relationship  Between  Pqq  and  pf  j 


In  equation  (11)  set  x  =  0  to  get, 

1  2+2k 


1=0  m=2 

(22) 

For  fc  =  0  we  get. 

(23) 

which  yields. 

i>f(o)  =  c^4m< 

(24) 

where 

_  1  ~  V^^o,o 

V2h^,i 

(25) 

Recall  from  above  that  the  p^k^  can  be  found  from  the  boundary  functions  evaluated 

at  the  boundary. 

Combine  this  information  with  equation  (24)  to  get. 

P\,i  (c  )  Po,o- 

(26) 

That  is,  we  only  need  to  evaluate  one  of  the  boundary  functions  at  the  boundary 
in  order  to  know  p^  ^  and  pf'j.  The  following  subsection  will  evaluate  the  boundary 
function  (f>Q{x)  at  x  =  0. 

3.2.3  Evaluation  of  the  Boundary  Function  <f>o{x)  at  x  =  0. 

In  this  subsection  the  straightforward  evaluation  of  the  boundary  function  ^o(x)  at 
X  =  0  will  be  given. 
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In  equation  (11)  let  A:  =  1  and  i  =  2  to  get, 

«if(2)  =  +  Af,4<^(0)).  (27) 

The  values  of  the  scaling  function  <^(z)  at  the  above  integers  are  <f>(0)  =  1/2(1  +  \/3) 
and  <^(1)  =  1/2(1  —  \/3)5  see  [9].  Using  these  numbers  and  the  coefficients  provided  in 
[3]  we  get  <^f(2)  =  —.3654971039.  Similar  to  above,  let  Ar  =  0  and  X  =  1  in  equation 
(11)  to  get 

=  v^/i^>f'(2)  +  v/2/ii2^(0),  (28) 

and  evaluate  to  get  <Po{i)  =  —1.1265992786.  In  this  way,  one  can  get  the  following 
values  for  the  boundary  functions:  <^f(l)  =  1.4058555636,  <?io(l/2)  =  .4123639530. 
Now  we  have  two  values  for  4>o{x)  at  i  =  1/2  and  x  =  1.  Furthermore,  the  two 
boundary  functions  are  lines  between  x  =  0  and  x  =  1.  We  can,  therefore,  find 
the  value  of  <f>o{x)  at  x  =  0:  <l>o{0)  =  1.9513271847.  As  outlined  in  the  previous 
subsection  we  now  have  sufficient  information  to  find  Pq  q  and 

Pop  =  -1.9038388908, 

and 

pf,  =  -.04295212263. 

The  calculation  of  the  remaining  elements  of  the  derivative  projection  matrix  is 
straightforward  and  will  not  be  given.  The  next  subsection  provides  the  complete 
matrix  D. 

3.3  The  Complete  Matrix  D 

This  section  of  the  paper  will  be  concluded  by  giving  a  10  by  10  version  of  the  matrix 
D.  Only  four  digit  accuracy  will  be  given  for  the  corner  entries  p,  a,  and  13.  The 
entries  r  are  known  exactly. 


D  = 


12 


Note  that,  as  expected,  away  from  the  boundaries  D  is  the  same  as  the  derivative 


projection  matrix  which  was  found  for  periodic  boundary  conditions,  see  [6].  Also, 
note  D  is  not  antisymmetric,  but  is  banded.  D  is  not  in  itself  the  differentiation 
matrix,  but  D  is  the  component  of  the  differentiation  matrix,  D  =  C~^DC,  which 
does  the  differencing.  Later  in  this  paper  it  will  become  obvious  how  one  can  simply 
‘look’  at  D  in  order  to  tell  the  maximum  differentiation  accuracy  at  the  boundaries. 
The  next  section  begins  the  discussion  of  choosing  a  quadrature  matrix,  C. 
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4  Approximating  Scaling  Function  Coefficients 


Once  again,  let  us  work  with  only  the  D4  wavelet.  Recall,  as  constructed  in  [3], 
there  are  two  boundary  functions  at  each  end  of  the  interval,  as  well  as,  the  usual  D4 
scaling  functions  away  from  the  boundaries.  Let  N  denote  the  number  of  degrees-oL 
freedom  of  the  approximation  space,  Vq.  Our  function  f{x)  €  L^{R)  is  then  discretized 
as  follows  (In  order  to  simplify  the  presentation  1  will  let  the  ‘interval’  be  [0,  A^], 
and  there  will  be  no  overlap  between  left-hand  boundary  functions  and  right-hand 
boundary  functions.) 

First,  the  left-hand  side, 

(29) 

Jo 

then  the  right-hand  side, 

=  /  <t>f{x)f{x)dx,  (30) 

Jo 

where  i  =  0, 1.  Next,  discretize  in  the  middle  with  the  usual  scaling  functions, 

Sk=  j  (f>{x  -  k)f{x)dx,  (31) 

Jo 

for  k  =  2, ...,  iV  —  4. 

Now  we  will  assume  that  f{x)  is  defined  on  a  discrete  evenly-spaced  grid,  /(x,) 
for  i  =  1,...,A^.  We  must,  therefore,  approximate  the  above  three  integrals  by  an 
appropriate  quadrature  matrix. 

In  the  following  two  sections  two  types  of  quadrature  matrices  will  be  consid¬ 
ered.  The  first  type  is  constructed  using  the  moments  of  the  scaling  function  and  is 
represented  by, 

s  =  Cf.  (32) 

It  will  be  seen  that  two-point  quadrature  matrices  for  this  method  are  well-conditioned, 
whereas  three-point  quadrature  matrices  are  ill-conditioned. 


The  second  type  of  quadrature  matrix  is  constructed  from  samples  of  the  scaling 
function  and  boundary  functions  and  is  represented  by, 

/=  Cs.  (33) 


It  will  be  seen  that  this  type  of  quadrature  matrix  is  always  ill-conditioned. 


15 


5  Quadrature  Matrix  C  Using  Moments 

The  quadrature  matrix  C  maps  from  evenly-spaced  point  values,  /(x^),  of  a  function 
to  the  approximate  scaling  function  coefficients  at  the  finest  scale,  i*.  The  quadrature 
matrix  depends  on  three  quantities:  the  number  of  vanishing  moments  of  the  wavelet, 
the  grid  that  the  function  f{x)  is  defined  on,  and  which  grid  points  are  used  to 
approximate  each  scaling  function  coefficient, 

5.1  Moments  of  Scaling  and  Boundary  Functions 

The  moments  of  the  scaling  functions  were  calculated  in  appendix  A  of  [6],  The 
first  3  moments  are,  Mq  =  1,  Mi  =  .6339745962,  and  M2  =  .4019237886. 

We  will  concentrate  here  on  calculating  the  moments  of  the  boundary  functions. 

5.1.1  Moment  0  of  Left-Hand  Side  Boundary  Functions 

The  0-th  moment  of  the  boundary  function  k  is, 

r<x>  /— .JL,  r  t  t 

4>l{x)dx  =  /  <^f(2x)  +  y/2J2  -  m).  (34) 

Now,  adopt  the  notation  mj  =  f  d>j^(x)dx  and  let  y  =  2x  where  appropriate  in  the 
above  integrals  to  get, 

1  2+2*  .00 

=  X]  +  X  I  (35) 

1=0  m=2 

The  integral  <l>{x)dx  is  1  since  m  >  2,  and  the  scaling  function,  as  defined  here, 

is  supported  on  [—1,2].  Equation  (35)  defines  the  following  2  by  2  system: 

^0,0  ~  Aq,!  rn®  _  ~^o,2 

^f,o  ^1,1  ~  .  .  ~^l,2  ~  ^1,3  ~ 

The  solution  to  the  above  system  is 

=  .3620520589, 
m°  =  1.001445402. 

In  a  similar  fashion  all  moments  can  be  found  for  the  left  and  right  hand  side 
boundary  functions. 
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5.2  The  Choice  of  Grid  Points 

In  this  subsection  the  affect  of  the  choice  of  different  grid  points  on  the  quadrature 
matrix  will  be  illustrated.  For  simplicity,  in  this  subsection  we  will  work  with  only 
8  ba^es  functions  and  8  grid  points.  The  interval  will  be  [0,8].  When  choosing  our 
grid  points  to  approximate  the  scaling  function  coefficient  s,  we  must  use  only  the 
grid  points  which  lie  within  the  support  of  the  basis  function  6,(x).  That  is,  the  inner 
product  that  is  being  approximated  is, 

=  /  b,{x)f{x)dx, 

Jo 

and  the  data  that  is  available  to  approximate  must  be  limited  to  function  values 
defined  within  the  support  of  6,(x).  The  support  of  each  of  the  eight  bases  functions 
is, 

supp{<f>Y}  =  [0,2], 
supp{</)f }  =  [0,3], 
suppih)  =  [1,4], 
supp{<i>4}  =  [2,5], 

5upp{<^5}  =  [3,6], 
supp{<l>e}  =  [4,7], 
supp{<f>^}  =  [5,8], 
supp{<f>l^}  =  [6,8], 

where  the  superscript  ‘fed’  denotes  a  boundary  function. 

The  following  subsections  will  consider  a  few  of  the  possible  quadrature  matrices 
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corresponding  to  the  following  two  grids: 


and 


.5  ■ 

1.5 

2.5 

3.5 


6.5 

7.5 


0 

8/7 

2*8/7 

3*8/7 

4*8/7 

5*8/7 

6*8/7 

8 


(37) 


(38) 


5.2.1  The  Entries  of  the  Quadrature  Matrices 

Now,  one  must  use  the  scaling  function  and  boundary  function  moments  and  at  least 
two  grid  points  contained  within  the  support  of  each  beisis  function  in  order  to  build 
the  quadrature  matrix.  The  first  row  of  the  quadrature  matrix  will  estimate  the  basis 
function  coefficient  for  the  boundary  function  with  the  smallest  support,  Suppose 
we  have  chosen  grid  x.  We  know  that  for  the  wavelet  that  constants  and  lines  can 
be  represented  exactly  in  the  wavelet  subspace  Vq.  There  are  only  two  grid  points 
within  the  support  of  Therefore,  the  solution  to  the  following  linear  system  will 
yield  the  two  coefficients  ci,i  and  ci,2  in  the  first  row  of  the  quadrature  matrix.  The 
linear  system  is. 


■  1  1  ■ 

Cl,l 

fj  (f)\^{x)xdx 

.5  1.5 

.  ^1-2  . 

In  a  similar  manner  all  the  entries  of  a  quadrature  matrix  can  be  found  once  one 
heis  calculated  the  moments  of  the  basis  functions  and  chosen  the  grid  points.  The 
following  subsection  will  illustrate  a  few  of  the  quadrature  matrices  associated  with 
various  grid  choices. 
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5.2.2  All  Grid  Points  Contained  in  Scaling  function  Support 


As  a  first  possible  choice  of  grid  points,  let  us  use  every  grid  point  that  is  contained 
within  the  support  of  the  respective  basis  function.  Such  a  quadrature  matrix  will 


the  following  structure: 


Cl,l 

Cl.2 

0 

0 

0 

0 

0 

0 

C2,I 

C2,2 

C2,3 

0 

0 

0 

0 

0 

0 

C3,2 

C3,3 

^3,4 

0 

0 

0 

0 

0 

0 

C4,3 

^4,4 

C4.5 

0 

0 

0 

0 

0 

0 

C5,4 

^5,5 

Cs.e 

0 

0 

0 

0 

0 

0 

C6,5 

C6,6 

C6,7 

0 

0 

0 

0 

0 

0 

C7,6 

C7,7 

C7,8 

0 

0 

0 

0 

0 

0 

C8,7 

C8,8 

(40) 


For  both  of  the  grids  x  and  y  the  matrix  is  ill-conditioned.  The  condition  number 
is  proportional  to  where  N  is  the  number  of  grid  points,  and  inversion  of 


unacceptably  corrupts  the  differentiation  matrix. 


5.2.3  2-Point  Quadrature  at  Boundary  and  3-Point  Otherwise 

The  quadrature  matrix  has  the  following  structure: 


^3pt  _ 


Ci.i 

Cl  .2 

0 

0 

0 

0 

0 

0 

0 

C2,2 

C2,3 

0 

0 

0 

0 

0 

0 

C3,2 

^3,3 

<^3,4 

0 

0 

0 

0 

0 

0 

C4,3 

<^4,4 

^4,5 

0 

0 

0 

0 

0 

0 

C5,4 

C5,5 

C5,6 

0 

0 

0 

0 

0 

0 

C6,5 

C6,6 

Cej 

0 

0 

0 

0 

0 

0 

C7,6 

CjJ 

0 

0 

0 

0 

0 

0 

0 

C8,7 

C8,8 

(41) 


This  matrix  is,  also,  ill-conditioned  for  both  grids  x  and  y.  Again,  the  condition 
number  if  proportional  to  N^. 
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5.2.4  A  First  Option  Using  a  2-Point  Quadrature 


= 


C],l  Ci_2  0 

0  C2,2  ^2,3 

0  0 

0 
0 
0 
0 
0 


0 

0 

0 

0 

0 


0  0 

0  0 

C3,3  ^34  0 

0  C4_4  C4S 

0  0 
0 
0 
0 


0 

0 

0 


0 
0 
0 
0 

C5.5  <^5,6 
0 
0 


0 

0 

0 

0 

0 


C6,6  Ce,7 
C7,6  C7  7 


0 

0 

0 

0 

0 

0 

0 


0  0  Cg  7  rgg 

This  quadrature  matrix  is  well-conditioned  for  both  grids  x  and  y. 


Second  Option  for  2-Point  Quadrature 

^1,: 

<^1,2 

0 

0 

0 

0 

0 

0 

C2,l 

<^2,2 

0 

0 

0 

0 

0 

0 

0 

0 

C3,3 

^3,4 

0 

0 

0 

0 

Cf  ‘  = 

0 

0 

0 

^4,4 

^4,5 

0 

0 

0 

0 

0 

0 

0 

C5.5 

^5,6 

0 

0 

0 

0 

0 

0 

0 

C6,6 

C6,7 

0 

0 

0 

0 

0 

0 

0 

C7,7 

C7,8 

0 

0 

0 

0 

0 

0 

Csj 

C8,8 

This  quadrature  matrix  is,  also,  well-conditioned  for  grids  x  and  y. 

5.3  Conclusion 


(42) 


(43) 


In  this  section  a  variety  of  structures  for  quadrature  matrices  have  been  shown  where 
the  matrices  are  found  by  using  the  moments  of  the  scaling  function.  In  general,  the 
pattern  is  that  3-point  quadrature  matrices  are  ill-conditioned  and  2-point  quadrature 
matrices  are  well-conditioned. 

The  next  section  will  explore  quadrature  matrices  constructed  from  the  samples 
of  the  scaling  functions. 
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6  Quadrature  Matrix  C  Using  Samples 


Let  g{x)  denote  the  projection  of  /(x)  6  in  Vq: 

gix)  =  PvJix)  =  YisM^)-  (44) 

fc=i 

The  quadrature  formula  based  on  samples  of  the  basis  functions  can  be  found 


g{x,)  =  Pvofixi)  =  Yi  ^kh{xi).  (45) 

fc=i 

In  matrix  form  the  quadrature  formula  is, 


g  =  Cl  (46) 


or  let  be  the  approximation  to  s  to  get, 

/  =  Ca,  (47) 


where  for  an  8  by  8  caise  the  quadrature  matrix  is, 

•  <^^(xi)  <^f'(xi)  0  0  0  0  0  0  ■ 

<f>o{^2)  <t>2{^2)  0  0  0  0  0 

0  <f>2(x3)  Mxa)  0  0  0  0 

0  0  <^2(x4)  <f>3{x4)  <f>4{X4)  0  0  0 

C  =  0  0  0  (i)3{xi)  4>4{x$)  <j>s{xi)  0  0  (48) 

0  0  0  0  (I>4{X6)  <f>b{xB)  <i>\{xe,)  0 

0  0  0  0  0  <^5(^7)  <f>\{x7)  (i>o[x7) 

0  0  0  0  0  0  <^f(x8)  <^?(X8) 


Now  C  will  be  found  for  the  two  previously  defined  grids  x  and  y. 


6.1  A  Quadrature  Matrix  based  on  Grid  x 

Finding  the  explicit  form  of  the  above  defined  C  for  the  grid  x  is  relatively  simple  since 
we  only  need  the  numerical  values  of  the  scaling  functions  and  boundary  functions 
halfway  between  the  integers. 
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6.1.1 

The  Numerical  Values  of  the  Scaling  Function 

Using  the  definition  of  the  scaling  function  supported  on  [0,3], 

3 

<f){x)  =  hk<t>{2x  -  k), 

k=Q 

we  get 

<^(5/2)  =  \/2h^<t>{2)  =  .06698729811, 

(49) 

<^(3/2)  =  x/2(/i,<^(2)  +  /i2<^(l)), 

(50) 

and 

<^(1/2)  =  \/2M(l)  =  .9330127019. 

(51) 

The  numerical  values  of  (f>  at  the  integers  are  provided  in  [9]: 

=  1/2(1 +n/3), 
and 

4>{2)  =  1/2(1  -  \/3), 

6.1.2  The  Numerical  Values  of  the  LHS  Boundary  Functions 

Given  the  numerical  values  of  the  scaling  function  ^{x)  on  grid  x  one  can  now  find 
the  numerical  values  of  the  boundary  functions  on  grid  x  from  the  expression, 

1  2+2fc 

+  V2  5:  hi„<K2x  -  m).  (52) 

1=0  m=2 

In  a  similar  manner  one  can  find  the  numerical  values  for  the  RHS  boundary 
functions  to  get  the  quadrature  matrix  of  the  next  subsection. 

6.1.3  The  Complete  Quadrature  Matrix  for  Grid  x 

An  8  by  8  example  of  the  quadrature  matrix  for  grid  x  is, 
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■  .4124  .8495  0  0  0  0  0  O' 

.2062  -.0077  .9330  0  0  0  0  0 

0  .0669  0  .9330  0000 

0  0  .0670  0  .9330  000 

Cj  =  0  0  0  .0670  0  .9330  0  0  .  (53) 

0  0  0  0  .0670  0  .8561  0 

0  0  0  0  0  .0670  .3270  .4451 

0  0  0  0  0  0  -.1406  .8902 


This  matrix  is  ill-conditioned  just  as  the  3-point  quadrature  formulas  were  for  the 
quadrature  matrix  derived  using  moments. 

6.2  A  Quadrature  Matrix  based  on  Grid  y 

Using  the  definition  of  the  scaling  function, 

3 

4>{x)  =  \/2  hk4>{2x  -  fc), 

Jk=0 

one  can  set  up  a  matrix  such  that  the  eigenvector  for  the  eigenvalue  A  =  1  gives  a  scalar 
multiple  of  the  scaling  function  <i>{x)  at  ar  =  r/7  for  i  =  1,...,20.  The  eigenspace  for 
the  eigenvalue  A  =  1  is,  however,  2-dimensional.  To  avoid  this  complication  one  can 
approximate  the  numerical  values  of  <f>{x)  at  i  =  i/7  for  i  =  1,...,20  by  generating 
<i>{x)  at  a  large  number  of  diadic  points,  in  this  case  2’^  points  were  chosen.  This 
approximation  is  sufficient  to  understand  the  character  of  the  quadrature  matrix, 
specifically  if  it  is  well-conditioned  or  not.  Once  the  numerical  values  of  the  scaling 
function  are  known  then  we  can  find  the  numerical  values  of  the  boundary  functions. 
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The  numerical  values  in  terms  of  the  scaling  function  values  on  grid  y  are, 


4{8n)  =  ^rKh^A{\^n)^Ka4>{:2ii)) 

<?if(8/7)  =  v^(/if,i^if'(l6/7)  +  Af.2<|i(2/7)  +  h'(^^4>{-T^ll))  (54) 

<^('(16/7)  =  ^{/ifXll/7)  +  /if'X4/7)). 

As  with  grid  x  the  quadrature  matrix  produced  using  grid  jf  is  ill-conditioned.  Re¬ 
call  from  the  previous  section  that  all  .3- point  quadrature  matrices  were  ill-conditioned. 
I,  therefore,  conclude  this  section  with  the  observation,  not  a  proof,  that  quadrature 
matrices  for  wavelet  bases  on  an  interval  using  the  construction  in  [3]  should  be  con¬ 
structed  using  moments  and  that  the  number  of  grid  points  used  to  approximate  each 
scaling  function  coefficient  should  be  equal  to  the  number  of  vanishing  moments  of 
the  corresponding  wavelet. 
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7  Order  of  Accuracy 

In  this  section  the  order  of  accuracy  for  a  few  of  the  differentiation  matrices  which  can 
be  constructed  from  the  derivative  projection  matrix  D  and  the  quadrature  matrices 
constructed  in  the  previous  three  sections  will  be  found.  We  begin  with  the  differ¬ 
entiation  matrices  constructed  using  grid  x  followed  by  the  differentiation  matrices 
constructed  using  grid  y. 

7.1  Accuracy  Using  Grid  x 

Recall  that  grid  x  is  comprised  of  the  points  halfway  between  the  integers  and,  con¬ 
sequently,  the  spacing  between  the  grid  points  is  the  same  as  the  translation  distance 
for  the  scaling  functions  on  the  finest  scale.  Furthermore,  this  grid  does  not  provide  a 
grid  point  at  the  boundary  itself.  This  lack  of  a  grid  point  at  the  boundary  makes  this 
grid  an  unlikely  choice,  but  I  have  chosen  to  study  the  accuracy  as  an  experiment. 
Let  us  begin  by  giving  a  few  explicit  examples  of  differentiation  matrices. 

7.1.1  Differentiation  Matrix  Grid  x 

Note  the  unusual  form  of  the  following  matrix.  It  is  essentially  upper  triangular  in 
form.  This  form  is  due  to  the  inversion  of  From  Table  (7.1.3)  it  can  be  seen  that 
this  differentiation  matrix  does  not  maintain  the  superconvergence  at  the  boundaries. 
That  is,  the  differentiation  is  Ist-order  accurate  at  the  boundaries,  wheretis  it  is  4th- 
order  accurate  at  the  middle  50%  of  the  grid  points.  The  total  differentiation  accuracy 
across  the  entire  domain  is,  therefore,  2nd-order. 

V  = 
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-25.33 

42.38 

-20.77 

3.73 

— 6e~'^ 

0 

0 

-7e-^ 

7e-3 

-9e-3 

2e-5 

36-3 

-5.71 

-1.59 

8.35 

-1.04 

-le-3 

0 

0 

-2e-^ 

.02 

-.02 

3e-® 

6e-3 

1.05 

-8.11 

.05 

8 

-1 

0 

0 

-3e-3 

.03 

-.04 

7e-5 

.01 

0 

1 

-8 

0 

8 

-1 

0 

-7e-3 

.07 

-.09 

26-“ 

.03 

0 

0 

1 

-8 

0 

8 

-1 

-.02 

.16 

-.20 

36-“ 

.06 

0 

0 

0 

1 

-8 

0 

8 

-1.03 

-.34 

-.43 

76-“ 

.13 

0 

0 

0 

0 

1 

-8 

0 

7.93 

-.28 

-.94 

2e-3 

.29 

0 

0 

0 

0 

0 

1 

-8 

-.16 

9.56 

-3.02 

36-3 

.62 

0 

0 

0 

0 

0 

0 

1 

-8.35 

3.36 

3.65 

-.99 

1.33 

0 

0 

0 

0 

0 

0 

0 

.25 

-.76 

-8.37 

6.02 

2.86 

0 

0 

0 

0 

0 

0 

0 

-1.62 

16.61 

-33.66 

11.96 

6.71 

0 

0 

0 

0 

0 

0 

0 

.53 

-5.47 

2.78 

-3.27 

5.43 

7.1.2  Differentiation  Matrix  ^  Grid  x 

Now  we  use  the  quadrature  matrix  But,  the  results  are  similar  to  the  previ¬ 

ous  matrix.  The  structure  is  unacceptable  and  the  superconvergence  is  lost  at  the 
boundaries  as  can  be  seen  from  Table  (7.1.3) 


V  = 


-12.21 

7.44 

10.48 

-6.44 

.74 

0 

0 

0 

0 

0 

0 

0 

-10.99 

-11.15 

44.98 

-24.53 

1.70 

0 

0 

0 

0 

0 

0 

0 

-2.50 

-1.01 

-3.50 

8 

-1 

0 

0 

0 

0 

.01 

-.01 

0 

.46 

.07 

-7.54 

0 

8 

-1 

0 

0 

0 

.01 

-.02 

.01 

0 

0 

1 

-8 

0 

8 

-1 

0 

-.03 

.02 

-.04 

.02 

0 

0 

0 

1 

-8 

0 

8 

-1 

-.01 

.05 

-.09 

.04 

0 

0 

0 

0 

1 

-8 

0 

8 

-1.01 

.11 

-.19 

.09 

0 

0 

0 

0 

0 

1 

-8 

0 

7.97 

-.75 

-.41 

.19 

0 

0 

0 

0 

0 

0 

1 

-8 

-.06 

8.53 

-1.89 

.04 

0 

0 

0 

0 

0 

0 

0 

1 

-8.13 

1.14 

6.09 

-.11 

0 

0 

0 

0 

0 

0 

0 

0 

.73 

-5.53 

-3.12 

7.93 

0 

0 

0 

0 

0 

0 

0 

0 

-.24 

2.74 

-16.75 

14.26 

7.1.3  Conclusion  for  Grid  x 


Recall  that  in  the  construction  of  the  differentiation  matrix  that  the  middle  matrix  D 
is  fixed.  The  Table  (7.1.3)  contains  the  accuracy  for  two  well-conditioned  quadrature 
matrices  and  In  addition,  as  an  experiment  1  have  combined  two  ill- 

conditioned  matrices:  one  mapping  from  f  to  s  and  the  other  mapping  from  s  to 
/• 
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Differentiation 

Grid 

Boundary 

Error 

Inner 

Error 

Total 

Error 

Matrix 

Size 

l\  Error 

Ratio 

l\  Error 

Ratio 

li  Error 

Ratio 

16 

.0066 

2.86e-“ 

.0024 

im 

Grid  X 

32 

.0032 

2.06 

3.20e-® 

89 

5.73e-“ 

64 

.0016 

2.00 

2.10e-® 

1524 

1.41e-“ 

BtlM 

128 

7.856-“ 

2.04 

2.80e-" 

75 

3.50e-^ 

4.03 

256 

3.81e-“ 

2.06 

1.75e-*2 

16.0 

8.53e-® 

4.10 

512 

1.69e-“ 

2.25 

1.07e-*3 

16.4 

1.94e-® 

4.40 

16 

.0091 

2.71e-® 

.0026 

Grid  X 

32 

.0047 

1.94 

3.13e-^ 

86.58 

6.78e-“ 

3.83 

64 

.0024 

1.96 

5.92e-i° 

529 

1.72e-“ 

3.94 

128 

.0012 

2.00 

2.80e-” 

21.1 

4.33e-® 

3.97 

256 

6.08e-“ 

1.97 

1.75e-’2 

16.0 

l.lOe-® 

3.94 

512 

3.22e-“ 

1.89 

1.07e-’^ 

16.4 

2.89e-® 

3.81 

CfDC?'" 

16 

.0038 

6.52e-5 

.0013 

Grid  X 

32 

.0019 

2.00 

1.60e-5 

4.08 

3.15e-“ 

4.13 

64 

9.23e-“ 

2.06 

3.95e-® 

4.05 

7.87e-5 

3.99 

128 

4.64e-“ 

1.99 

9.83e-'^ 

4.02 

1.97e-® 

3.99 

256 

2.32e-“ 

2.00 

2.45e-^ 

4.01 

4.92e-^ 

4.00 

512 

1.18e-“ 

1.97 

6.12e-« 

4.00 

1.22e-® 

4.04 

Table  1;  Differentiation  accuracy  for  Daubechies  wavelets  defined  on  an  interval  using 
grid  points  halfway  between  the  integers. 


The  noteworthy  points  are,  first  of  all,  that  in  all  three  cases  the  accuracy  at  the 
boundary  is  first  order.  Recall  that  when  periodic  boundary  conditions  are  imposed 
that  the  differentiation  accuracy  is  4th  order  for  the  D4  wavelet.  This  loss  of  su¬ 
perconvergence  is  a  serious  problem.  Furthermore,  note  from  the  table  that  the  two 
well-conditioned  quadrature  matrices  maintain  the  superconvergence  away  from  the 
boundary.  But,  in  all  three  cases  the  total  error  of  the  differentiation  matrix  across 
the  entire  domain  is  one  higher  than  the  boundary  accuracy,  i.e.,  we  get  2nd  order 
differentiation  accuracy. 

Note:  Inner  accuracy  is  calculated  at  the  middle  50%  of  the  grid  points,  and 
boundary  accuracy  is  calculated  at  the  two  outermost  grid  points  at  each  end  of  the 
interval. 
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7.2  Accuracy  Using  Grid  y 

We  now  study  the  differentiation  accuracy  where  we  have  included  a  boundary  point 
at  each  end  of  the  interval  in  the  grid  and  the  remainder  of  the  grid  points  are 
required  to  be  evenly-spaced  across  the  interval.  This  choice  of  grid  is  a  much  more 
likely  choice  than  grid  x.  However,  superconvergence  is,  again,  lost  at  the  boundary. 
This  reduces  the  total  accuracy  of  the  differentiation  to  2nd  order  as  with  grid  x. 
Again,  let  us  note  a  few  explicit  examples  of  differentiation  matrices. 

7.2.1  Differentiation  Matrix  Grid  y 

Again  note  the  unusual,  and  undesirerable,  form  of  the  following  matrix.  As  can  be 
seen  from  Table  (7.2.2),  again,  the  superconvergence  is  lost  at  the  boundary. 


D  = 


-1.99 

3.33 

-1.60 

.26 

-26-“ 

-le-“ 

-66-3 

1 

1 

l6-“ 

-26-“ 

66-3 

26-3 

-.46 

-.14 

.72 

-.13 

-le-3 

-66-“ 

-36-“ 

1 

1 

76-“ 

-16-3 

36-“ 

l6-“ 

.08 

-.63 

-.02 

.70 

-.12 

-26-3 

-16-3 

-7e-“ 

2e-3 

-46-3 

16-3 

36-“ 

0 

.07 

-.63 

-.02 

.69 

-.12 

-36-3 

-26-3 

76-3 

-.01 

36-3 

16-3 

0 

0 

.08 

-.63 

-.02 

.69 

-.11 

-56-3 

.02 

-.03 

96-3 

36-3 

0 

0 

0 

.08 

-.63 

-.02 

.69 

-.12 

.05 

-.08 

.02 

76-3 

0 

0 

0 

0 

.08 

-.63 

-.02 

.68 

5e-3 

-.18 

.05 

.02 

0 

0 

0 

0 

0 

.08 

-.63 

-.04 

.94 

-.50 

.12 

.03 

0 

0 

0 

0 

0 

0 

.08 

-.69 

.52 

-.12 

.15 

.07 

0 

0 

0 

0 

0 

0 

0 

-.03 

.43 

-1.54 

1.00 

.14 

0 

0 

0 

0 

0 

0 

0 

-.20 

2.12 

-4.25 

2.05 

.29 

0 

0 

0 

0 

0 

0 

0 

.44 

-4.62 

7.00 

-2.83 

86-3 

It  has  been  mentioned  many  times  that  the  structure  is  undesirable.  This  is  due 
to  the  lack  of  directional  independence  of  the  matrix.  One  ‘fix’  for  this  problem  is 
given  in  the  next  subsection. 

7.2.2  The  NPD  Differentiation  Matrix  Grid  y 

NPD  denotes  No-Preferred-Direction.  It  is  possible  to  fix  one  of  the  problems  with 
the  above  matrices,  the  problem  of  directional  dependence.  That  is,  a  general  differ¬ 
entiation  matrix  should  never  treat  data  coming  from  the  right  differently  than  data 
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coming  from  the  left.  This  problem  can  be  ‘fixed'  if  one  simply  turns  the  basis  around 
and  builds  another  differentiation  matrix.  That  is,  turn  all  the  inner  scaling  functions 
around  and  interchange  the  role  of  the  RHS  and  LHS  boundary  functions.  The  new 
differentiation  matrix  T)  can  be  found  from  the  unreversed  differentiation  matrix  V 
by  first  flipping  the  matrix  by  a  middle  horizontal  axis  and  then  by  a  middle  vertical 
axis  and  then  negating.  That  is, 

i)  =  -flipud{fliplr{V)),  (55) 


and  we  build  a  new  NPD  differentiation  matrix  by  averaging  these  two  matrices: 

'Dnpd  =  \{T>  +  V).  (56) 

The  following  NPD  matrix  treats  data  coming  from  the  left  in  the  same  manner 
as  data  coming  from  the  right.  Note  that  this  matrix  is  full. 


-12.00 

36.96 

-51.58 

29.30 

-2.67 

-7e-* 

-3e~* 

-2e-* 

8e-* 

-le-3 

4e-=' 

le-" 

-4.47 

-13.12 

29.85 

-13.47 

1.22 

-4e-» 

-2e-’ 

-le-3 

4e-3 

-7e-® 

2e-3 

6e-* 

-.39 

-9.77 

9.11 

1.62 

-.54 

-.01 

-6e-® 

-4e~® 

.01 

-.02 

7j-3 

2e-3 

-.43 

-.43 

-3.01 

-3.24 

8.30 

-1.16 

-.02 

-.01 

.04 

-.07 

.02 

6e-3 

-.20 

-.70 

3.44 

-9.44 

.12 

7.96 

-1.15 

-.03 

.11 

-.18 

.06 

.02 

-.09 

-.32 

1.09 

.42 

-7.86 

-9e-® 

7.94 

-1.16 

.28 

-.47 

.14 

.04 

-.04 

-.14 

.47 

-.28 

1.16 

-7.94 

9e-® 

7.86 

-.42 

-1.09 

.32 

.09 

-.02 

-.06 

.19 

-.11 

.03 

1.15 

-7,96 

-.12 

9.44 

-3.44 

.70 

.20 

-6e-^ 

-.02 

.07 

-.04 

.01 

.02 

1.16 

-8.30 

3.24 

3.01 

.43 

.43 

-2e-3 

-7e-3 

.02 

-.01 

4e~® 

6e-3 

.01 

.54 

-1.62 

-9.11 

9.77 

.39 

— 6e~^ 

-2e-3 

7e-3 

-4e-3 

le-3 

2e-® 

4^—3 

-1.22 

13.47 

-29.85 

13.11 

4.47 

-le-* 

— 4e“'^ 

le-3 

-Se-* 

2e-* 

36-* 

7e-* 

2.67 

-29.29 

51.58 

-36.96 

12.00 

Table  (7.2.2)  gives  the  accuracy  results  for  two  differentiation  matrices  constructed 
from  the  well-conditioned  quadrature  matrix 

Again,  the  following  table  shows  that  loss  of  superconvergence  at  the  boundaries 
regardless  of  the  differentiation  matrix.  Also,  note  that  for  grid  y  that  the  supercon¬ 
vergence  is  not  even  achieved  away  from  the  boundary  as  it  is  for  grid  x.  Inner  error 
and  boundary  error  are  defined  in  the  same  manner  as  for  the  previous  table. 

In  conclusion,  in  this  section  I  have  presented  a  number  of  differentiation  matrices 
along  with  the  orders  of  accuracy  of  the  differentiation.  In  all  cases  the  superconver¬ 
gence  which  is  encountered  under  periodic  boundary  conditions  is  lost. 
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Differentiation 

Grid 

Boundary 

Error 

Inner 

Error 

Total 

Error 

Matrix 

Size 

1\  Error 

Ratio 

l\  Error 

Ratio 

/i  Error 

Ratio 

16 

.0181 

8.70e-^ 

.00598 

Grid  y 

32 

1.94 

2.316-“ 

3.77 

.00162 

3.69 

64 

6.04e-5 

3.82 

4.26e-“ 

3.80 

128 

.00238 

1.99 

1.49e-® 

4.05 

1.089e-“ 

3.91 

256 

.00117 

2.03 

3.73e-^ 

3.99 

2.72e-® 

4.00 

NPD  of 

16 

6.82e-“ 

.00692 

32 

1.42e-® 

48.0 

.00191 

3.62 

Grid  y 

64 

.00606 

1.91 

5.88e-« 

2.41 

5.05e-“ 

3.78 

128 

.00310 

1.95 

1.43e-® 

4.11 

1.30e-“ 

3.88 

256 

.00156 

1.99 

3.50e-’' 

4.09 

3.29e-* 

3.95 

Table  2;  Differentiation  accuracy  for  Daubechies  wavelets  defined  on  an  interval  us¬ 
ing  grid  points  at  the  boundaries  and  all  other  grid  points  evenly-spaced  across  the 
interval. 
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8  Other  Constructions  and  Conclusion 


In  this  paper  I  have  explored  the  differentiation  matrix  for  Daubechies-type  wavelets 
defined  on  an  interval  using  the  boundary  construction  defined  in  [3].  The  exploration 
has  been  limited  to  the  wavelet  with  2  vanishing  moments  usually  referred  to  as 
the  Daubechies  4  wavelet.  It  wa.s  seen  that  the  superconvergence  encountered  with 
periodic  boundary  conditions  which  was  proven  in  [6]  is  lost  with  this  boundary 
construction. 

There  are  other  boundary  constructions  available  for  Daubechies  wavelets,  see  [1], 
[8],  but  none  of  the  constructions  can  maintain  the  superconvergence.  There  is  a 
very  straightforward  way  to  find  the  maximum  accuracy  at  the  boundary  for  a  given 
boundary  construction.  That  is,  when  the  basis  functions  are  orthogonal  under  trans¬ 
lation,  as  they  are  for  Daubechies  wavelets,  one  can  simply  count  the  number  of  basis 
functions  which  have  support  overlapping  the  support  of  the  outermost  boundary 
functions.  This  defines  tlie  number  of  nonzero  entries  in  the  first  and  last  rows  of  the 
middle  matrix  D  referred  to  here  as  the  derivative  projection  matrix.  All  differencing 
in  wavelet  differentiation  is  done  by  D,  and,  therefore,  the  differentiation  accuracy 
at  the  boundaries  cannot  exceed  one  less  than  the  number  of  nonzero  entries  in  the 
first  and  last  rows  of  D.  Of  course,  this  argument  holds  for  every  row  of  D,  but  I  am 
particularly  concerned  with  the  boundaries  here. 

Perhaps  it  is  not  possible  to  construct  boundary  functions  for  wavelets  which  give 
everything.  That  is,  boundary  functions  which  maintain  the  wavelet  structure,  which 
maintain  the  approximation  properties  across  the  entire  interval  and  which  maintain 
superconvergence.  This  conjecture  is  based  first  on  the  fact  that  to  date  there  is  not 
a  boundary  construction  which  satisfies  these  criteria  for  either  Daubechies  wavelets 
or  spline-based  wavelets,  see  [7].  In  further  support  of  this  conjecture,  in  a  paper 
by  Gottlieb,  Gustafsson,  Olsson,  and  Strand,  [5],  it  is  proven  that  it  is  not  possible 
to  construct  boundary  functions  for  a  finite  element  subspace  which  maintain  the 


superconvergence  encountered  with  periodic  boundary  conditions  if  there  are  charac¬ 
teristics  leaving  the  domain.  Note  that  the  finite  element  subspace  is  the  same  as  the 
first  order  spline  wavelet  subspace,  see  [7]. 

In  conclusion,  under  the  currently  available  boundary  constructions  the  supercon¬ 
vergence,  which  is  encountered  under  periodic  boundary  conditions,  is  not  maintained 
at  the  boundaries.  In  the  case  of  the  D4  wavelet  under  the  boundary  construction 
examined  here  one  gets  4th-order  differentiation  away  from  the  boundaries  and  only 
Ist-order  differentiation  at  the  boundaries  producing  a  scheme  with  only  2nd-order 
accurate  spatial  differentiation.  It  remains  to  be  seen  if  it  is  possible  to  construct 
boundary  functions  for  a  wavelet  subspace  which  can  maintain  superconvergence  at 
the  boundaries. 
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